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Abstract. We study the structural and thermodynamic properties of a model of 
point particles interacting by means of a Gaussian pair potential first introduced by 
Stillinger [StilHnger F H 1976 J. Chem. Phys. 65 3968]. By employing integral 
equation theories for the fluid state and comparing with Monte Carlo simulation results, 
we establish the limits of applicability of various common closures and examine the 
dependence of the correlation functions of the liquid on the density and temperature. 
We employ a simple, mean-field theory for the high density domain of the liquid and 
demonstrate that at infinite density the mean-field theory is exact and that the system 
reduces to an 'infinite density ideal gas', where all correlations vanish and where the 
hypernetted chain (HNC) closure becomes exact. By employing an Einstein model 
for the solid phases, we subsequently calculate quantitatively the phase diagram of 
the model and find that the system possesses two solid phases, face centered cubic and 
body centered cubic, and also displays reentrant melting into a liquid at high densities. 
Moreover, the system remains fiuid at all densities when the temperature exceeds 1% 
of the strength of the interactions. 

Published: J. Phys.: Condens. Matter 12, 5087 (2000) 
PACS numbers: 61.20.-p, 61.20.Gy, 64.70.-p 

1. Introduction 

The structural and phase behaviour of systems whose constituent particles interact 
by means of pair potentials diverging at zero separation between them, is a problem 
that has been extensively studied in the last three decades, by a variety of theoretical, 
experimental and computational methods. It has by now been established that the 
diverging repulsions at close separations between the particles are the dominant factor 
causing crystallisation of the system; indeed, for the purposes of understanding the 
freezing mechanism, the simple, hard sphere model is in most cases sufficient. The 
crystal structure into which a system freezes depends on the steepness of the repulsion, 
with hard repulsions favouring a face centered cubic (fee) lattice and soft ones a 
body centered cubic (bcc) lattice [^. On the other hand, interparticle attractions are 
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responsible for bringing about a liquid-gas coexistence, whose stability with respect to 
the freezing transition depends crucially on their relative range with respect to that of 
their repulsive counterparts 0. 

Though repulsions are necessary to bring about a solidification transition, they 
are by no means sufficient. This was demonstrated recently by Watzlawek et al ^ 
who studied a system interacting by means of an ultrasoft, logarithmically diverging 
potential which has been shown to model accurately the effective interaction between 
star polymers in a good solvent [^, |^. In particular, the pair potential employed in Refs. 
iSl, 01, reads as 



Mr) = ^/^/^ 



In ' 



a 



v7(r - cr) 
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for r < (J 



for r > cr, 



where /3 = {ksT)'^ is the inverse temperature {ks is Boltzmann's constant), / is the 
functionality of the stars and cr is the typical extension of the stars, the so-called corona 
diameter. It was found that the strength of the repulsion, controlled by the parameter / 
in eq. is crucial in determining whether the system crystallises. For values / < 34, 
the system remains fluid at all densities and, even for values / > 34, reentrant melting 
into a high-density liquid has been observed for a particular range of functionalities 

We are, therefore, faced with a question which is specular to the old question of how 
much attraction is needed in order to make a liquid 0, namely: How much repulsion is 
needed in order to make a solid? In this respect, it is interesting to consider the extreme 
family of potentials that are hounded, i.e., they remain finite for the whole range of 
interparticle separations, even at full overlap between the particles. In the context 
of microscopic interactions in atomic systems, such bounded potentials are evidently 
unphysical: full overlaps between atoms, molecules or even compact macromolecular 
entities are forbidden by the Born repulsions between the electrons. Yet, they are 
perfectly realistic in the context of effective interactions between fractal objects such as 
polymer chains, as will be demonstrated in section 0. 

In comparison with systems interacting by means of unbounded potentials, very 
little is known about bounded interactions. A model system that has been recently 
studied is the 'penetrable spheres model' (PSM), in which the interaction is a positive, 
finite constant for separations smaller than a diameter cr and vanishing otherwise 
The model was studied by means of cell-model calculations and computer simulations 
P], liquid-state integral equation theories and density- functional theory and 
it was found that no reentrant melting takes place with increasing density because a 
clustering mechanism stabilises the solid at all temperatures. The PSM is, however, 
rather peculiar in this respect and both the constant value of the potential inside cr 
and its sharp fall to zero outside cr render it into a rather unrealistic model for effective 
interactions. 
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A much more realistic model of a bounded potential is the Gaussian core model 
(GCM), introduced in the mid-seventies by Stillinger [jlT|. In the GCM, the interaction 
potential between the particles reads as 

t;(r) = £e-('-/'^)', (2) 

where £ > is an energy and a a length scale. Stillinger's original work was 



complemented later by molecular dynamics simulations fl^, |13|, high-temperature 



expansions and the discovery of exact duality relations in the crystalline state [15 



Based on these considerations, a semi-quantitative phase diagram of the GCM was 



drawn |Tl], 0, which displayed reentrant melting and a maximum freezing temperature 
t^. The latter is defined as the temperature above which the system remains fluid at 
all densities. Yet, a detailed study of the structure of the GCM fluid by means of the 
modern techniques of liquid-state theory and a quantitative calculation of the phase 
diagram of the GCM are still lacking. The purpose of this work is to fill this gap and to 
draw from the study of the GCM some general conclusions, applicable to a large class 
of bounded potentials. 

The rest of the paper is organised as follows: in section ^ we briefly review previous 
work on this model and make contact with effective interactions between polymer chains. 
In section ^ we present results about the structure of the fluid, obtained by means of 
liquid-state integral equation theories and Monte Carlo simulations. In section |^ we 
discuss the high density limit of the fluid, by employing a simple mean-field density 
functional. In section |^ we present the theory for the crystal phases and in section ^ we 
make use of the results for the fluid and the solids in order to quantitatively draw the 
phase diagram of the system. Finally, in section ^ we summarise and conclude. 

2. Gaussian effective interactions 

A Gaussian pair potential was first proposed fifty years ago by Flory and Krigbaum 



17|, as the effective interaction between the centres of mass of two polymer chains. 



This so-called Flory-Krigbaum potential, vfk{i^), reads as 

where Vseg and fsoiv denote the volumes of a monomer segment and a solvent molecule, 
respectively, is the degree of polymerisation, Rg the radius of gyration of the chains 
and X is a parameter that controls the solvent quality (0 < x < 1/2 denotes a good 
solvent and x > 1/2 a poor one.) 

Formally, the effective interaction potential fefr('") between two polymer chains is 
defined as follows: the centres of mass of the two chains are held fixed at separation 
r, and the canonical trace is taken over all the monomer degrees of freedom. In this 
way, a restricted partition function Q{r) is obtained, where the restriction denotes the 
constraint of holding the centres of mass fixed. The effective interaction is then simply 
the free energy —kBTlnQ{r). As can be seen in Fig. 0, the centres of mass of two 
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polymer chains can coincide without any violation of the excluded volume condition 
between the monomers (beads). Hence, the effective interaction between polymer centres 
of mass does not have to diverge at r = and a bounded potential is fully realistic. 



Figure 1. A snapshot from a simulation involving two self-avoiding polymers. In 
this configuration, the centres of mass of the two chains (denoted by the big sphere) 
coincide, without violation of the excluded-volume conditions. (Courtesy of Arben 
Jusufi.) 

The effective interaction between linear polymers has been extensively studied by 



computer simulations, both on- lattice |18, ff9l and off-lattice [EG, B^, E2l E3l 24 1 and 



for a varying solvent quality. In addition, Kriiger et al ||2^ performed a theoretical 
calculation of the effective interaction for athermal solvents, using renormalisation-group 
techniques. Regardless of the detailed manner in which the monomer-monomer excluded 
volume interactions have been modeled, all these studies converge to the result that 
the effective potential fcfr(r) is extremely well approximated by a universal Gaussian 
function^ given by eq. (0), with e being an energy scale of the order of the thermal 
energy ksT and a a length scale of the order of the spatial extent of the chains (e.g., 
the radius of gyration). This type of effective Gaussian interaction has been recently 
employed by Louis et al ||25| in their study of the structure of colloid-polymer mixtures. 

The Gaussian effective potential between polymer chains in purely entropic in 
nature, causing the energy e to scale with temperature. Therefore, the latter completely 
drops out of the problem as a relevant thermodynamic parameter. However, since the 
system is of significant theoretical interest, it is pertinent to consider e and fc^T in more 
generality, as two independent energy scales and study the properties of the system 
as functions of density and temperature. To this goal, we define the dimensionless 
temperature t as 

t^^-^ = m-\ (4) 

The density p of a system of particles (or polymer chains in the above discussion) 
enclosed in the macroscopic volume V is the ratio N/V. As dimensionless measures of 
the density, we are going to use to use in this paper both the parameter p = pcr^ and 
the 'packing fraction' rj, defined as 

V = ^P(^^- (5) 
o 

A number of exact properties for the GCM were worked out by Stillinger et al 
TT| , |15|, |l^ , yielding important information on the topology of the phase diagram. The 



relevant properties are listed below. 

% The universality holds when the two chains have the same degree of polymerisation and the solvent is 
athermal. If the temperature becomes relevant or the two chains have different numbers of monomers, 
corrections to the Gaussian form appear. 
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1. Hard sphere limit. At low temperatures and densities, the GCM reduces to a hard 
sphere (HS) system with a HS diameter that diverges as the temperature approaches 
zero [jlT|. Hence, the GCM displays there the usual HS freezing transition from a fluid 
into a fee lattice and the coexistence densities converge to zero at vanishing temperature. 
Using the known results for the HS freezing p6| , p^ , p8| the shape of the freezing and 



melting lines tf{p) and tm(p) as t ^ is given by the equations 

(6) 
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Up) = i^exp(^--^J. (7) 

2. Duality relations. Due to the property of the Gaussian potential to remain form- 
invariant under Fourier transformation, it has been shown that the internal energy of a 
certain crystal at low densities (lattice sum) is related to that of the reciprocal crystal 



at high densities ||T5[. In this way it was established that the fee and bcc lattices are 
degenerate at t = at the 'magic density' p^: = Tr~^^'^ = 0.1796, with the fee "winning" 
for p < p<t: and the bcc for p > p^. 

3. Reentrant melting. At small but finite temperatures and at high densities the 
bcc lattice remelts into a fluid and the high density freezing and melting lines are given 
by the equations |T5[ 

tf(p) ocexp(-irfp2/3), (8) 
tm(p) oc exp (-iri^p2/3) , (9) 

with the appropriate constants and K^. Note the duality between eqs. (|]), (|^), valid 
for p — > 0, and (|), ®, valid for p — > cxd. Apart from proportionality constants, one 
pair can be obtained from the other by the formal substitution p 1/p. Moreover, 
eqs. (H) and (^) show that the slope of the high-density freezing and melting lines are 
negative. Contrary to the usual case, the liquid coexisting with the solid has a higher 
density than the latter. 

The above general considerations were combined with molecular dynamics 
simulations [|l^, [13| and an approximate phase diagram of the model was drawn [l^ . 



It was found that below an 'upper freezing temperature' t^ ~ 0.008, the GCM displays 
the transitions fluid fee, fee bcc and bcc fluid with increasing density, whereas 
above tu a single fluid phase exists and no freezing takes place, at any density. However, 
the liquid structure of the GCM has not been studied to-date and the phase diagram 
was drawn only semi-quantitatively. In the following sections we perform a quantitative 
analysis on all these aspects, using standard tools of statistical mechanics. 

3. The fluid: integral equation theories and Monte Carlo simulations 

The theoretical determination of the pair structure of a uniform fluid amounts to the 
calculation of the radial distribution function g{r) and the direct correlation function 



(dcf) c(r) p9|. The correlation function h{r) is simply g{r) — 1 and then h{r) and c(r) 
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are connected to each other through the Ornstein-Zernike relation which has the form 
2| 



h{r) = c(r) ~^ P J dr'c{\r — v'\)h{r') 



(10) 



With h{Q) denoting the Fourier transform of h{r), the structure factor S{Q) is defined 
as 

S{Q) = l + ph{Q). (11) 

The Ornstein-Zernike relation is exact. As it connects two unknown functions, one 
more relation or closure is needed in order to determine g{r) and c(r). Closures are 
approximate relations which arise from exact diagrammatic expansions of g{r) in terms 
of c(r) but with certain classes of diagrams ignored. The exact relation between g{r) 
and c(r) reads as 



g{r) = exp [-pv{r) + h{r) - c{r) + 5(r)] , (12) 

where v{r) is the pair potential and B{r) is the so-called bridge function, consisting of 
the sum of all elementary diagrams that are not nodal. 

All known closures can be thought of as approximate relations for the form of 
B{r). Common closures are the Percus-Yevick (PY) and Hypernetted Chain (HNC) 
approximations In the PY closure, the approximation for B{r) reads as 

5pY(r) = c(r) - h{r) + ln[^(r) - c(r)], (13) 

whereas in the HNC the approximation is made that the bridge function vanishes: 

5HNc(r) = 0. (14) 

The PY closure is successful only for short-ranged, hard interactions and hence it 
will not be considered here. In addition to the HNC, we have considered two more 



closures of increasing degree of sophistication, the Rogers- Young (RY) closure and 
the zero-separation (ZSEP) closure ||31], 

as 

1' 



g{r) = exp [—(3v{r)] 



1 + 



, |3^, |3^. The former (RY) closure reads 
exp [-f{r)f{r)] 



fir) 



where 



7(r) = h{r) — c(r) 
and the 'mixing function' /(r) is chosen to have the form 
/(r) = 1 — exp (— ar) . 



(15) 



(16) 



:i7) 



The parameter a is determined so that thermodynamic consistency between the virial 
and compressibility pressures is achieved In the limit a = one recovers the PY 
and in the opposite limit, a — > oo, the HNC closure. 
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The ZSEP closure includes three parameters and is a direct approximation for the 
bridge function having the form 



(18) 



1 + Q;7(r) 

The ZSEP closure has been applied recently by Fernaud et al |^ to the PSM model 
mentioned above, yielding excellent results for the correlation functions as compared 
to the simulation results. The three parameters '-P and a are determined in such a 
way that virial-compressibility, Gibbs-Duhem and zero separation consistencies are all 
fulfilled; for details on the above we refer the reader to Ref. 0. Evidently, for C = the 
ZSEP reduces to the HNC. 

In addition, in order to test the reliability of the various closures, we have performed 
standard NVT Monte Carlo (MC) simulations for different densities and temperatures. 
In what follows, we focus our attention to temperature t = 0.01 on two grounds: 
on the one hand, according to the approximate phase diagram |6[, at t = 0.01 the 
system remains fluid at all densities and therefore we can study the development of the 
correlation functions with density for an unlimited range of the latter, without entering 
a region where the fluid is metastable. On the other hand, this temperature is low 
enough, so that significant structure is expected for the correlation functions of the fluid 
and hence the integral equation theories can be put in a strong test. We present the 
obtained results in Figs. ||, and ^ and we discuss them below. 
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Figure 2. Comparison for the radial distribution function g{r) between the simulation 
results and those obtained from the various closures, at t = 0.01 and at small values 
of the packing fraction, (a) 77 = 0.05; (b) 77 — 0.12. 



We begin with the RY closure. When the packing is not large enough, typically 
77 < 0.50 at this temperature, the RY closure gives results which are in very good 
agreement with simulations, as can be seen in Figs. Q and |^. However, above 77 ~ 0.50, 
the g{r) from simulations starts penetrating towards r = 0, physically meaning that the 
probability of finding two particles 'sitting on top of each other' is finite and there is no 
'correlation hole' in g{r). This is to be expected for a system with a bounded interaction. 
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However, as can be seen from Fig. ^ (a), the RY closure fails to capture precisely this 
penetration phenomenon, yielding g{r)^s that are too low at small separation and making 
the erroneous prediction that a correlation hole exists. 

The reason for this behaviour can be traced back in the construction of the RY 
closure, eqs. (|T5]) and (p!?]), where it can be seen that the RY closure always looks like 
the PY closure at small separations. The latter is however inaccurate for a long-range 
interaction lacking a hard core. We attempted to modify the RY closure by employing 
mixing functions yielding a HNC-like small-r behaviour and a PY-like large-r behaviour. 
However, this did not yield self-consistency parameters for the whole range of densities. 
Despite of its inability to correctly describe the high-density penetration of g{r), the 
standard RY closure yields nevertheless a self-consistency parameter a for all densities 
considered here. In addition, this parameter grows with density, thus pointing to a 
tendency of the RY closure to reduce to the HNC at high packings. 





r/a 



Figure 3. Same as Fig. g but for intermediate packings, (a) rj = 0.20; (b) rj — 0.50. 
Note the different scale of the vertical axes. 




Figure 4. Same as Fig. ||but for high packings, (a) rj = 1.00; (b) 77 — 6.00, where the 
simulation result is indistinguishable from the HNC. 
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In view of the failure of the RY-closure, we are led to consider the ZSEP closure 
which has precisely the property that in its implementation the value of g{r) at zero 
separation plays an important role and is determined self-consistently. In fact, the 
resounding success of the ZSEP in describing the g{r) of the PSM model (also a bounded 
interaction) has been mainly attributed to this property 0. As can be seen from Figs. 
^ and ^(a), the ZSEP performs only slightly worse than the RY closure up to a packing 
fraction 7] ^ 0.25. The self-consistency parameters (, (p and a of the ZSEP are displayed 
in Fig. H. The parameter C, which appears as an overall factor in the parametrisation of 
the bridge function [see eq. ([18|) ] decreases with density and close to ?7 = 0.25 it is small 
enough and the ZSEP practically reduces to the HNC closure. 

However, at packing fractions rj > 0.25, a second branch of solutions appears, which 
is denoted by the dashed lines in Fig. ^. This branch is disconnected from the first and 
hence it causes the g{r) to behave discontinuously at this packing fraction, a result 
which is clearly unphysical. This second branch produces (y'(r)'s that show too much 
penetration and too little structure, when compared with the simulation results. The 
reason for this unphysical behaviour can be traced to the fact that this second branch 
corresponds to bridge functions which are positive at small separations r. Indeed, from 
eq. (0), and taking into account that 7(r) ^ 1 at small separations, we can see that 
-BzsepI'^) has the sign of the product ({ip — a) /a. For the first, physical branch (sohd 
lines in Fig. |^) this combination is negative because ip < a. For the second, unphysical 
branch, this combination is positive because ip > a there. A positive bridge function acts 
then as an additional 'attractive potential' in eq. ( ]T^ ) and causes the overpenetration 
in g{r) mentioned above. 




0.25 0.5 0.75 1 1.25 1.5 1.75 2 



Figure 5. The self-consistency parameters of tlie ZSEP closure applied to the Gaussian 
core model at reduced temperature t — 0.01 as functions of the packing fraction rj. The 
solid lines denote the physical and the dashed lines the unphysical branch (see the text.) 
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The appearance of a branch of solutions of the ZSEP for which the bridge function 
is positive is a sign of internal inconsistency of the closure and in this sense the ZSEP 
signals its own limits of applicability. Indeed, the bridge function of any system has 
been shown to be a quasi-universal function, which can always be mapped onto the 
bridge function of a suitably defined hard sphere system having a hard sphere diameter 
that depends on the characteristics of the interaction and the thermodynamic point 
under consideration |^. The bridge function of the HS system is, however, practically 
exactly known and it is essentially a negative function for all r. Hence, a positive bridge 
function is physically unacceptable and the second branch of solutions of the ZSEP 
has to be discarded. The results from this second branch come again into very good 
agreement with simulation at packing fractions above r] ^ 1.00 because, as seen in Fig. 
^, the parameter ( is already very small there and the bridge function has a negligible 
contribution to g{r); even the unphysical branch reduces to the HNC at high densities. 
However, there is no way to bridge the physical solutions at packings t] < 0.25 with the 
HNC-like solutions within the ZSEP closure, that is without having to compare with 
independently produced simulation results. 

We comment next on the quality of the HNC. As can be seen from Figs. - 
the HNC underestimates the structure at small to intermediate packings but yields 
otherwise reasonable results. It does not suffer from any of the problems of the more 
sophisticated closures and, in fact, it seems to be the best at high densities. In order to 
further explore this property (which is supported by the fact that the other two closures 
tend to the HNC at this limit), we also solved this closure at extremely high packing 
fractions. Here, the highest packing at which we simulated was rj = 6.00, due to time 
constraints. With increasing rj, a very large number of particles would be required in 
the simulation box in order to obtain reliable results. In view of the fact that the HNC 
gives results which coincide with those from Monte Carlo at r/ = 6.00 [see Fig. §(b)], 
we gain confidence at this closure and apply it for arbitrarily high densities, in order to 
obtain information on the structure of the very dense fluid. In Fig. ^(a), we show results 
for g{r) where it can be seen that at very high packings g{r) tends to unity and hence 
h{r) = g{r) — 1 tends to zero. However, this does not mean that the structure factor 
S{Q) tends to unity as well, as a naive guess based on the definition S{Q) = 1 + ph{Q) 
would imply. The quantity h{Q) tends to zero but at the same time the prefactor p 
diverges, thus giving rise to a nontrivial S{Q). 

Results for the corresponding structure factor S{Q) are shown In Fig. ^b). It can 
be seen that, for high densities, the peak of S{Q) disappears and the latter looks like a 
'smoothed step function' with values ranging from zero to one. The value of Q at which 
the crossover occurs does not scale as a power-law of the density but rather shifts to the 
right by almost a constant every time the packing fraction is increased by an order of 
magnitude. This hints to a very weak dependence of this crossover value on the density; 
we return to this point in section ^. 

The liquid-state correlation functions of the GCM display an anomalous behaviour 
in comparison with that of 'normal' liquids, interacting by means of hard, unbounded 




Figure 6. (a) The radial distribution function g{r) and (b) the structure factor S{Q) 
as obtained by the HNC closure for extremely high values of the packing fraction. 



interactions. For the latter, the structure grows monotonically with increasing density 
and eventually the systems freeze. Here, the structure grows up to a packing fraction 
7] = 0.12 at t = 0.01 and for higher densities it becomes weaker again. To amply 
demonstrate this phenomenon we show in Fig. ^ the structure factors obtained from the 
MC simulations for a wide range of densities, where it can be seen that the height of 
the peak of S{Q) attains its maximum value at rj = 0.12. The same phenomenon has 
been observed in the liquid structure of star polymer solutions employing the interaction 
potential given by eq. (P p, ||, ^ and, in fact, in star polymer solutions the effect has 
also been observed experimentally |33]. This behaviour of S{Q) is closely related to 
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Figure 7. The structure factor S{Q) of the GCM at t 
fractions, as obtained by Monte Carlo simulation. 
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reentrant melting 0. Indeed, the height of the maximum of S{Q) is a diagnostic tool 
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for the freezing transition. According to the Hansen- Verlet criterion a hquid 



crystalhses when S{Q) at its highest peak has the quasi- universal value 2.85. At higher 
values of the peak the system is solid and at lower fluid. The evolution of S{Q) with 
density shown in Fig. |^ in conjunction with the Hansen- Verlet criterion implies that the 
system approaches crystallization at about rj = 0.12 and then remelts. The height of 
the peak at ?7 = 0.12, which is indeed slightly above the Hansen- Verlet value, implies 
that at this temperature the system barely freezes and that t = 0.01 is very close to the 
upper freezing temperature. We will confirm this prediction in section ^ where we will 
also use the structural results obtained here in order to calculate the free energy of the 
fluid. 

Finally, we have performed MC simulations and solved the HNC closure also at 
much higher temperatures, t = 1.00, corresponding to the physical domain for the 
effective interactions between polymer chains. There, we found that the HNC accurately 
reproduces the simulation results at all densities and that the liquid has very little 
structure, a result which can be easily understood in view of the fact that the thermal 
energy, which is equal to the interaction strength there, washes out the correlation effects 
caused by the interactions. 

4. The high-density limit in the fluid state 

The results of the preceding section on the structure of the fluid at high densities point 
out that the correlations are becoming weaker as density grows and that the system 
approaches some kind of 'infinite density ideal gas' limit, where g{r) = 1. This limit 
was assumed already by Stillinger There, the internal energy of a high-density fluid 
was approximated by that of a random distribution of points interacting by means of 
the Gaussian potential, i.e., the positions of the points were assumed to be uncorrelated. 
The relation g{r) = 1 was used explicitly there in deriving an estimate for the internal 
energy of the high-density fluid and it was shown that in fact the complete absence 
of a correlation hole in the fluid raises its internal energy with respect to that of a 
solid (where a correlation hole is necessarily present) by exactly e/2 per particle. Yet, 
the existence of this ideal-gas limit was not proven. Here, we are going to present a 
density functional mean-field theory which establishes this limit and provides analytic 
expressions for the correlation functions of the GCM fluid at high densities. 

In the framework of density functional theory (DFT), one considers in general a 
system with a spatially modulated one-particle density p(r). The Helmholtz free energy 
of the system is then a unique functional of the one-particle density , F = F[p{r)], 



and can be written as a sum of an entropic, ideal part Fid[p(r)] which is exactly known 
and an excess part Fex[p(r)] which is in general unknown: 

F[p(r)]=Fid[p(r)]+Fe.[p(r)] 

= ksT J dvp{r) [ln[p(r)A3] - l] + Fejp(r)]. (19) 



In eq. (|l^) above, A denotes the thermal de Broglie wavelength of the particles. 
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We consider the hmit pa^ ^ 1. The average interparticle distance a = p~^/^ 
becomes vanishingly small in this limit and it holds a <^ a, i.e., the potential is extremely 
long-range. Every particle is simultaneously interacting with an enormous number of 
neighbouring molecules and in the absence of short-range excluded volume interactions 
the excess free energy of the system can be approximated by a simple mean-field term, 
equal to the internal energy of the system: 



^ex[p(r)] = IJ J drdr'v{\r-r'\)p{r)p{r') 



(20) 



with the approximation becoming more accurate with increasing density and eventually 
exact at p ^ oo. The direct correlation function c(|r — r'|; p) in a fluid of density p is 
given by the second functional derivative of -Fex[p(r)] with respect to the density ||41|| , 
namely 

c(r-r ;p) = - lim (21) 
P(r)-P 6p{r)6p{r') 



Combining eqs. ( pOD and ( pTD we find that the dcf of the system at high densities becomes 
independent of p and is simply proportional to the interaction: 

c(r) = -pv{r). (22) 

Eq. ( P^ ) above has a strong similarity with the mean-spherical approximation 
(MSA) introduced as a perturbation theory to study systems interacting by 

potentials that can be separated into a hard sphere interaction (diameter o"hs) ^o('") and 
a 'soft tail' 0(r). In the MSA, the radial distribution function of the system g{r) vanishes 
for r < (Ths and the direct correlation function c(r) at r > cths is given by c(r) = — /?0(r). 
The main difference from the MSA here, is that there is no reference potential fo(r) 
because there are no hard cores in the system. Hence, eq. (^) holds with the total 
interaction on the right hand side and for the whole range of separations r. Moreover, 
unlike the MSA which is essentially a high-temperature approximation, eq. (|2^ ) holds 
for the whole temperature range, provided that the density is high enough {pa^ ^1). Of 
course, the mean- field approximation becomes also valid at high temperatures {t ^ 1) 
irrespective of the density, because there the thermal energy completely dominates over 
the bounded interaction. This happens in contradistinction with diverging interactions, 
where short-range correlation effects always survive, at all temperatures. 

The Fourier transform c{Q) of the dcf is obtained easily from eqs. ( [2^ ) and (0), and 
for the GCM it has the form: 

d{Q) = -7c^/^f3£a^ exp [- {QcT/2f] . (23) 

Using the Ornstein-Zernike equation and the ensuing relation S{Q) = [1 — pc(Q)]^^ we 
obtain for the structure factor of the GCM the expression: 

S(Q) = — ^. (24) 

This analytic expression is compared with the MC result at (3e = 100 {t = 0.01) and 
T] = 6.00 in Fig. |^. The excellent agreement between the two demonstrates the validity 




Figure 8. The structure factor S{Q) of the GCM att = 0.01 and rj = 6.00 as obtained 
from simulation (soHd Une) and as given by the analytical expression, eq. (^) (dashed 
line). 



of the simple mean-field theory at high densities. Eq. (^) shows that at high densities 
S{Q) is a monotonia function of Q and has the shape of a 'smoothed step function' with 
values ranging from (vr^/^/Jepcr^) = at low Q's to unity at high Q's. The crossover 
between the two regimes occurs at a characteristic wavenumber at which S{Q) = 1/2 
and which, according to eq. (pi), is given by 



Q,a = 2^1n(7r3/2/?£pa3). (25) 

Note the very weak, square root-logarithmic dependence of Q* on density and the inverse 
temperature. 

Another quantity of interest is the isothermal compressibility xt of the system, 
defined as 

^r-{v^,) , (26) 

where F is the Helmholtz free energy, and also related to the Q ^ limit of S{Q) 
through 

pkBTxT = S{0). (27) 



The expressions ([19|) and (^) yield the Helmholtz free energy in the high density limit 



as 

1 



F = ksTN [In (pA^) - l] + -Nn^/^pa^, (28) 



and from eqs. (EBf) and (53) we obtain 



_ 1 

~ ksTp + 7r3/%pV3 ' ^^^^ 
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which evidently satisfies the compressibihty sum rule ( p7D with S{Q) given by eq. (p^). 
We note that at high densities xt obeys the scaling 

ksTxT ~ tp-\ (30) 

Moreover, from eq. ( p8| ) the pressure P = —{dF/dV) is obtained as 

P = ksTp + ^vr^/^^pV^. (31) 

Since S{Q) = 1 + ph{Q), eq. (0) immediately yields an analytic expression for the 
Fourier transform h{Q) of the correlation function h{r), namely 

1 + TT'^/^pepcT'^ exp [— [Qa/2) \ 

At low Q's, where the exponential is of order unity, the term proportional to the density 
in the denominator dominates and h{Q) behaves as — p^^ ^ as p — oo. At high Q's, 
the exponential itself tends to zero. Hence, the function h{Q) vanishes as p ^ oo with 
the leading term being proportional to p^^. Consequently, the correlation function h{r) 
vanishes as well and g{r) ^ 1 as p — oo. This is the particular 'high density ideal gas' 
limit of the model. However, a clear distinction must be drawn between this 'interacting 
ideal gas' and the usual ideal gas, in which either the system is noninteracting or an 
interacting system is considered at the opposite limit, p — > 0. In the usual ideal gas 
limit, we have c(r) = exp[— /3t;(r)] — 1, S{Q) = 1 and g{r) = exp[— /?t>(r)]. The ideal gas 
pressure P scales linearly with the density and the ideal compressibility xt scales with 
p~^. Here, c(r) = —i3v{r), S{Q) is not unity in the whole Q range, the pressure scales as 
p^ [see eq. (0)] and the isothermal compressibility as p~^ [see eq. (|30D]. Nevertheless, the 
above considerations point out to a kind of interesting duality of the GCM in the liquid 
phase, in which the system is trivially ideal at low densities and becomes again ideal 
(vanishing correlations) at high densities. This can be thought of as the counterpart for 
the fiuid state of the duality discovered by Stillinger in the crystalline state ||T5|] . 

The mean-field approximation put forward in this section is not limited to the 
Gaussian potential. It should be valid for all interactions v{r) which are finite, 
analytic functions and which tend to zero fast enough at infinite separations so that 
the thermodynamic functions (e.g., the internal energy) are extensive and the Fourier 
transform of v{r) exists. Hence, we make the theoretical prediction that this peculiar 
high-density ideal gas limit should exist for all potentials satisfying the requirements put 
forward above. In addition, it can now be understood why the HNC showed the best 
agreement with simulation results among all closures at high densities. As g{r) — >■ 1, 
h{r) and c(r) — —j3v{r) in this limit, the exact relation ([T^) forces the bridge 
function B{r) to vanish. For bounded, analytic potentials decaying fast enough to zero, 
the HNC becomes exact at high densities. 
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5. The solid: Einstein model 

In this section we present the approach employed for the calculation of the free energies 
of candidate crystalline states of the model, necessary for determining the phase diagram 
of the GCM. As we are dealing with a soft interaction, a harmonic approximation in the 



solid is justified and we adopt the simple Einstein model ^ as a means to estimate 
the free energy of the latter. The approach is based on the Gibbs-Bogoliubov inequality 
||29| . The latter states that the Helmholtz free energy F of a system having Hamiltonian 
7i is related to the Helmholtz free energy Fq of a reference system having Hamiltonian 
Ho by 

F<Fo + {n- no)o, (33) 

where the canonical average on the right hand side is taken in the reference system. 
The procedure is useful if (i) a simple enough reference Hamiltonian can be chosen, 
which physically corresponds to a situation close enough to the real one, and in which 
Fq and the average {Ti — 7io)o can be calculated in a straightforward way and (ii) this 
Hamiltonian contains at least one variational parameter which can be chosen in order 
to minimise the right hand side, obtaining in this way a lower upper bound for the true 
free energy of the system. 

In a harmonic solid, a reasonable non-interacting reference system is the Einstein 
solid, characterised by the Hamiltonian 



N r 2 1 



+ - (rj - Rj 

2m 2 ^ ' 



(34) 



where p, is the canonical momentum of a particle of mass m, k is the 'spring constant' 
which plays the role of a variational parameter and the set {Rj} forms a prescribed 
Bravais lattice. For the GCM, the real Hamiltonian reads as 

«-E|t+EE-p -V^ ■ (35) 

i=l i=l j=i+l ^ ' 

The calculation of the Helmholtz free energy of the Einstein solid is a trivial exercise 
yielding 

^ = |fc^Tlnf?V3A:BTlnf^V (36) 

where 



N 2 Vvr/ \a 



a = . (37) 

The one- and two particle densities in the reference system are given by a sum of 
Gaussians and a double sum of products of Gaussians, respectively, where the latter are 
centered on lattice sites and the sums are carried upon all these sites. As the interaction 
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has itself a Gaussian form, the calculation of the quantity {H — TCo)o reduces to integrals 
that can be carried out analytically. The final result reads as 

The sum on the rhs is carried over all shells of lattice vectors, i.e., sets of lattice vectors 
of equal length, with the shell j = (self-interaction) excluded. The quantity Uj is the 
number of lattice vectors belonging to a shell and Xj = |Rj|/cr. It can be easily seen 
that at T = 0, where a — >■ cxo, the first term on the rhs of eq. ( pSj) reduces to the internal 
energy per particle of the considered crystalline arrangement, the lattice sum. At zero 
temperature there is no variational parameter and the Einstein model becomes exact; 
the winning phase is the one with the lowest lattice energy. At finite temperatures, the 
sum of the terms on the rhs of eqs. (|36| ) and (|38| ) has to be minimised with respect to 
a for any given lattice structure. The minimum comes about through the competition 
of the entropic, logarithmic term on the rhs of eq. (|36|), which favors delocalisation at 
(5 = and the internal energy term on the rhs of eq. ( p8| ) which favors localisation at 
a oo. The obtained value is then the estimate for the Helmholtz free energy of the 
given lattice and the procedure can be repeated for every candidate lattice. We note 
that the term 3fcsTln(A/(j) on the rhs of eq. (^) can be dropped because it occurs for 
all possible phases of the system, fluid and solid, and does not affect the free energy 
comparisons between them. 

We have performed the minimisation for a different number of candidate lattice 
structures: fee, bcc, simple cubic, diamond and body-centered orthogonal, in which the 
ratios between the lattice constants of the conventional unit cells were used as additional 
variational parameters |0, We always assumed a lattice with a single occupancy per 
site, based on the result of Stillinger stating that indeed solids with multiply occupancies 



are unstable |]TT|. For the whole range of temperatures < t < 0.015 considered, we 
always found the fee and bcc to be the only stable crystals, with the former winning 
at low densities and the latter at high densities. The quantitative results are shown in 
section |6|. 



6. Quantitative phase diagram 

In this section we switch from the variable t] to the variable p = pa^ as a measure 
of the density, in order to allow for a direct comparison with the approximate results 
of Stillinger and Stillinger |TB[. With the free energy of the crystals obtained by the 
procedure outlined in section |^, the phase diagram can be drawn if the corresponding 
free energy of the fluid is also known. The latter can be obtained from the results of 
the integral equation theories outlined in section ||. 

The Helmholtz free energy of the liquid is the sum of the ideal and excess terms 
(see section |), namely 

F{p,T) = NksT [In {pa') - l] + SArfc^Tln (-) + F,^{p,T). (39) 
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If for a fixed temperature the radial distribution function g{r) is known for a region of 
densities, then one possibihty to calculate the free energy is through the so-called virial 
route. Here, one calculates for every density the excess pressure Pex of the system via 

Pc. = ^ / r'v\r)gir)dr, (40) 



3 

where v'{r) = dv{r)/dr. Then, the excess free energy can be calculated by integrating 
the thermodynamic relation Pgx = —dFi,^/dV from p = up to the given density, under 
the initial condition -Fex(p = 0, T) = 0. Alternatively, one can use the results for the 
structure factor S{Q) and the definition of the isothermal compressibility xt [eqs. (^) 
and (|27D ] to calculate the liquid free energy, following in this way the compressibility 
route. 

If thermodynamic consistency between the two routes is not explicitly enforced 
in an otherwise approximate closure, then the results from two routes are different, a 
problem known as thermodynamic inconsistency of the closure [^. The RY- and ZSEP- 
closures are thermodynamically consistent but the HNC is not. However, as explained 
in section ^, neither the RY nor the ZSEP yield reliable results for the whole density 
range. If we were dealing with a usual system, without reentrant melting, then we 
could have used the RY or the ZSEP results at low or intermediate densities. However, 
we are interested also in the high-density fluid free energies, where in fact the HNC 
becomes exact. A combination of low-density results from one closure and high-density 
results from another is not of much use either, because it would produce unphysical 
discontinuities in the free energy or its derivatives at the point of switching between the 
two. We are thus led to employ the HNC closure in the whole density domain in order 
to perform the thermodynamic integration and obtain the fluid free energy. 

The HNC compressibility route yields at low densities fluid free energies that are 
too low, leading to the erroneous result that the HS-like freezing of the Gaussian fluid 
into a fee lattice does not take place. There are two factors playing an important role 
here: on the one hand, the predicted isothermal compressibilities are too high causing 
a fluid free energy which is too low and on the other hand, the solid free energy, being 
a product of the variational procedure outlined in section ^, is unavoidably higher than 
the true one, see eq. (|33D . It is therefore pertinent to follow the HNC virial route in 
calculating the fluid free energy. The latter leads indeed to an overestimation of the 
fluid free energy but this compensates for the overestimation of the solid free energy 
and leads to the physically correct picture of freezing into an fee solid at the low-density 
part of the phase diagram. We have thus calculated the fluid free energies through the 
HNC virial route for a range of temperatures 10~^ t < 0.015 for a range of densities 
< pcr^ < 1.00 and performed the common tangent construction on the resulting 
Ffiuid(p)/^- and Fsoiid(p)/V"-curves to obtain the phase boundaries. The resulting phase 
diagram is shown in Fig. |^. 

The phase diagram obtained is very similar to the approximate one drawn by 
Stillinger and Stillinger |l^. It shows the sequence of freezing, structural (fee — >■ bcc) 
and remelting transitions as well as the upper freezing temperature t^ associated with the 



Fluid and solid phases of the Gaussian core model 19 
0.012 



0.01 



0.008 



0.006 



0.004 
0.002 



"0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 

3 

pa 

Figure 9. The phase diagram of the GCM obtained by the approach described in the 
text. The fcc-bcc coexistence hnes are also double lines but they cannot be resolved 
in the scale of the figure because the fcc-bcc density gap is too small. The full dot 
marks the point at which the fluid-bcc coexistence curves turn around. The two insets 
show details of the phase diagram, (a) In the neighbourhood of zero densities and 
temperatures, (b) In the neighbourhood of the fluid-fcc-bcc triple temperature, with 
the dashed line denoting the triple line between these coexisting phases. 

corresponding density pu- The coordinates of this point, where the fluid-bcc coexistence 
hnes turn around, are (tu,Pu) = (0.0102,0.2292). This is in perfect agreement with the 
prehminary results from section |^, where at t = 0.01 the structure factor at rj = 0.12, 
corresponding to p = 0.2292, was found to slightly exceed the Hansen- Verlet value. The 
fcc-bcc coexistence lines run linearly from the points pfcc = 0.17941 and pbcc = 0.17977 
at t = to the points pfcc = 0.16631 and pbcc = 0.16667 at the triple temperature 
tt = 8.753 X 10~^. The density gap between the fee- and bcc-coexisting densities remains 
constant and equal to Ap = 3.6 x 10~^. The density of the coexisting fluid at the triple 
temperature is pfluid = 0.16431. 

It should be emphasised that, notwithstanding its deceiving appearance in Fig. |, 
the point (tu,Pu) is not a critical point [|ri|]. At (tu,pu), two common tangents between 
the fluid- and bcc-solid free energies, one lying on the low- and the other on the high- 
density side of it, coalesce into this single point. No susceptibility diverges and all free 
energy density curves remain strictly concave up. 

It is now pertinent to ask whether the Hansen- Verlet freezing criterion is satisfied 
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for both the low- and the high-density crystalhsation of the system. To this end, we 
have performed additional MC simulations at temperatures below t^^ and in Fig. |10| we 
show structure factors at two such temperatures, t = 0.007 and t = 0.005, for increasing 
values of the density. It can be seen that the Hansen- Verlet criterion is indeed valid for 
both freezing transitions, a feature also observed for the reentrant melting phenomenon 
in star polymer solutions P, ffl. 
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Figure 10. Structure factors S{Q) for the GCM at (a) t = 0.007 and (b) t = 0.005 
obtained from MC simulations. Densities as indicated in the legends. The straight 
lines mark the Hansen- Verlet value 2.85. The corresponding structure factors in the 
regions 0.153 < pcr^ < 0.363 for t = 0.007 and 0.114 < pa^ < 0.477 for t = 0.005 show 
Bragg peaks, indicating that the GCM has to be in the solid state for these densities. 



It is of course also possible to calculate the exact free energies of the various 
candidate phases in a simulation, by means, e.g., of the virial route in the fluid state and 
by employing the Frenkel-Ladd method [0, ^, P, ^] in the solid state. However, the 
latter is very time-consuming. The very good agreement between the phase boundaries 
obtained from the approximate theory presented here and the MC results regarding the 
height of the peak of S{Q) and the spontaneous crystallisation of the system in the 
simulation box, give us confidence that the phase diagram of Fig. ^ is quantitatively 
correct. 



7. Summary and concluding remarks 

The Gaussian core model displays a whole range of unusual phenomena and properties: 
an anomalous dependence of the correlation functions on the density, a high-density 
'ideal gas limit', reentrant melting and an upper freezing temperature. Many of these 
characteristics arise from the fact that the pairwise interaction does not diverge at zero 
separation, i.e., it is bounded. However, the topology of the phase diagram of the GCM 
has some striking similarities with that of star polymers [0, obtained by employing 
the diverging interaction given by eq. An upper freezing parameter (t in the GCM 
and / in the star potential) appears in both, above which the systems remain fluid 
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at all densities. Freezing and reentrant melting occur also at both systems. Yet, the 
phase diagram of the GCM displays only two solid phases, whereas that of stars has a 
richness of exotic crystal structures 0. The latter are caused by the crossover of the 
star potential from a logarithmic to a Yukawa form, a feature absent in the GCM. 

To the best of our knowledge, the only other bounded interaction that has been 
studied by similar techniques to-date is the penetrable spheres model (PSM) mentioned 
in section ^ P, ^ ITof. The phenomenology and the associated phase diagram of the 
PSM are quite different from those of the GCM. These two bounded potentials have 
phase diagrams which do not look at all similar to one another. This comes as a 
striking difference to the relative insensitivity of the phase diagrams of unbounded 
interactions. In the PSM no reentrant melting occurs and the system seems to freeze at 
a// temperatures into increasingly clustered solids P, |T^. The reason for this difference 
is, evidently, that in the PSM, the particles can build clusters. In the GCM, this 
mechanism is not present; the interaction varies rapidly enough with distance, so that 
multiple occupancies are penalised. Therefore, the boundedness of the interaction brings 
about a new factor to be considered: multiple occupancies, which are prohibited from 
the very beginning for diverging interactions, have to always been taken into account 
whenever one deals with bounded ones. 

An important conclusion which appears to be valid for a large class of bounded 
potentials has been nevertheless drawn, and it is the mean field-ideal gas behaviour of 
all such systems at high densities. The PSM does not belong to this class; indeed, the 
discontinuity of the PSM pair potential at r = cr forces its radial distribution function 
to have a jump at the same position at all densities and hence an ideal gas limit can 
never be attained in this model. This is a feature complementary but not identical 
to the clustering property of the system. On the other hand, one can easily construct 
bounded pair interactions depending on some parameter, so that the sharpness of the 
decay of the PSM from a finite to zero value can be tuned. It will be very interesting 
to examine the structure and thermodynamics of such a family of systems as a function 
of the 'smoothening parameter' and establish the limits of the clustering- and the ideal 
gas-behaviour at high densities. We plan to return to this problem in the future. 
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